#include "math_integral.h"
#include <stddef.h> // use size_t
#include <cassert>
#include <algorithm>
#include <functional>
#include <cmath>
#include "constants.h"

namespace ModuleBase
{

Integral::Integral(){}

Integral::~Integral(){}


// Peize Lin accelerate 2017-10-02
/*
void Integral::Simpson_Integral
(
    const int mesh,
    const double *func,
    const double *rab,
    double &asum
)
{
    //     simpson's rule integration. On input:
    //     mesh = mhe number of grid points (should be odd)
    //     func(i)= function to be integrated
    //     rab(i) = r(i) * dr(i)/di * di
    //     For the logarithmic grid not including r=0 :
    //     r(i) = r_0*exp((i-1)*dx) ==> rab(i)=r(i)*dx
    //     For the logarithmic grid including r=0 :
    //     r(i) = a(exp((i-1)*dx)-1) ==> rab(i)=(r(i)+a)*dx
    //     Output in asum = \sum_i c_i f(i)*rab(i) = \int_0^\infty f(r) dr
    //     where c_i are alternativaly 2/3, 4/3 except c_1 = c_mesh = 1/3
    
    //  simpson's rule integrator for function stored on the
    //  radial logarithmic mesh
    //	routine assumes that mesh is an odd number so run check
    if (mesh % 2 == 0)
    {
        std::cout << "\n error in subroutine simpson ";
        std::cout << "\n routine assumes mesh is odd but mesh = "
             << mesh << std::endl;
        return;
    }

    asum = 0.00;
    const double r12 = 1.00 / 12.00;
    double f3 = func [0] * rab [0] * r12;
    for (int i = 1;i < mesh;i += 2)
    {
        const double f1 = f3;
        const double f2 = func [i] * rab [i] * r12;
        f3 = func [i + 1] * rab [i + 1] * r12;
        asum += 4.00 * f1 + 16.00 * f2 + 4.00 * f3;
    }
    return;
}// end subroutine simpson
*/


// Peize Lin accelerate 2017-10-02
void Integral::Simpson_Integral
(
    const int mesh,
    const double * const func,
    const double * const rab,
    double &asum
)
{
    /*     simpson's rule integration. On input:
    !      mesh = mhe number of grid points (should be odd)
    !      func(i)= function to be integrated
    !      rab(i) = r(i) * dr(i)/di * di
    !      For the logarithmic grid not including r=0 :
    !      r(i) = r_0*exp((i-1)*dx) ==> rab(i)=r(i)*dx
    !      For the logarithmic grid including r=0 :
    !      r(i) = a(exp((i-1)*dx)-1) ==> rab(i)=(r(i)+a)*dx
    !      Output in asum = \sum_i c_i f(i)*rab(i) = \int_0^\infty f(r) dr
    !      where c_i are alternativaly 2/3, 4/3 except c_1 = c_mesh = 1/3
    */
    //  simpson's rule integrator for function stored on the
    //  radial logarithmic mesh
    //	routine assumes that mesh is an odd number so run check
    assert(mesh&1);

    asum = 0.00;
	const size_t end = mesh-2;
    for( size_t i=1; i!=end; i+=2 )
    {
		const double f1 = func[i]*rab[i];
		asum += f1 + f1 + func[i+1]*rab[i+1];
    }
	const double f1 = func[mesh-2]*rab[mesh-2];
	asum += f1+f1;
	asum += asum;
	asum += func[0]*rab[0] + func[mesh-1]*rab[mesh-1];
	asum /= 3.0;
    return;
}// end subroutine simpson


// Peize Lin accelerate 2017-10-02
void Integral::Simpson_Integral
(
    const int mesh,
    const double * const func,
    const double dr,
    double &asum
)
{
    /*     simpson's rule integration. On input:
    !      mesh = mhe number of grid points (should be odd)
    !      func(i)= function to be integrated
    !      rab(i) = r(i) * dr(i)/di * di
    !      For the logarithmic grid not including r=0 :
    !      r(i) = r_0*exp((i-1)*dx) ==> rab(i)=r(i)*dx
    !      For the logarithmic grid including r=0 :
    !      r(i) = a(exp((i-1)*dx)-1) ==> rab(i)=(r(i)+a)*dx
    !      Output in asum = \sum_i c_i f(i)*rab(i) = \int_0^\infty f(r) dr
    !      where c_i are alternativaly 2/3, 4/3 except c_1 = c_mesh = 1/3
    */
    //  simpson's rule integrator for function stored on the
    //  radial logarithmic mesh
    //	routine assumes that mesh is an odd number so run check
    assert(mesh&1);

    asum = 0.00;
	const size_t end = mesh-2;
    for(size_t i=1; i!=end; i+=2 )
    {
		const double f1 = func[i];
		asum += f1 + f1 + func[i+1];
    }
	const double f1 = func[mesh-2];
	asum += f1+f1;
	asum += asum;
	asum += func[0] + func[mesh-1];
	asum *= dr/3.0;
    return;
}// end subroutine simpson


// Peize Lin add 2016-02-14
void Integral::Simpson_Integral_0toall
(
    const int mesh,
    const double * const func,
    const double * const rab,
    double * const asum
)
{
    // asum(r) = \int_{r'=0}^{r} dr' f(r') 

    const double r2=1.00/2.00, r3=1.00/3.00;
    asum[0] = 0.00;
    double f3 = func [0] * rab [0];
    for( int i=1; i<mesh; i+=2)
    {
        const double f1 = f3;
        const double f2 = func[i] * rab[i] ;
        f3 = func[i+1] * rab[i+1] ;
        asum[i] = asum[i-1] + r2*( f1 + f2);
        if(i+1<mesh)
        {
            asum[i+1] = asum[i-1] + r3*( f1 + 4.00*f2 + f3 );
        }
    }
    return;
}


// Peize Lin add 2016-02-14
// faster but still have bug
/*void Integral::Simpson_Integral_alltoinf
(
    const int mesh,
    const double *func,
    const double *rab,
    double *asum
)
{
    // asum(r) = \int_{r'=r}^{+\infty} dr' f(r') 
    //         = \inf_{r'=r}^{mesh} dr' f(r')

    const double r2=1.00/2.00, r3=1.00/3.00;
    asum[mesh-1] = 0.00;
    const int odd_mesh = (mesh-1)^~1;
    double f1 = func[odd_mesh] * rab[odd_mesh];
    for( size_t i=(mesh-3)|1; i>0; i-=2)
    {
        const double f3 = f1;   
        if( i+3==mesh )
        {
            const double f4 = func[mesh-1] * rab[mesh-1];
            asum[mesh-2] = r2*(f3 + f4);
        }
        const double f2 = func[i] * rab[i] ;
        f1 = func[i-1] * rab[i-1] ;
        asum[i-1] = asum[i+1] + r3*( f1 + 4.00*f2 + f3 );
        asum[i] = asum[i-1] - r2*( f1 + f2);
    }
    return;
}*/


// Peize Lin add 2016-06-11
// a little lower
void Integral::Simpson_Integral_alltoinf
(
    const int mesh,
    const double * const func,
    const double * const rab,
    double * const asum
)
{
    Integral::Simpson_Integral_0toall( mesh, func, rab, asum );

    const double asum_all = asum[mesh-1];
    for (int i = 0;i < mesh; ++i)
	{
        asum[i] = asum_all - asum[i];
	}
	return;
}

double Integral::simpson(const int n, const double* const f, const double dx) 
{
    assert(n >= 2);
    if (n == 4)
    { // Simpson's 3/8 rule
        return 3.0 * dx / 8 * (f[0] + 3.0 * f[1] + 3.0 * f[2] + f[3]);
    }

    if (n == 2)
    {
        return 0.5 * dx * (f[0] + f[1]);
    }

    if (n % 2 == 1)
    { // composite Simpson's 1/3 rule
        double sum = 0.0;
        for (int i = 1; i != n-2; i += 2) 
        {
            sum += 2.0 * f[i] + f[i+1];
        }
        sum += 2.0 * f[n-2];
        sum *= 2.0;
        sum += f[0] + f[n-1];
        return sum * dx / 3.0;
    }
    else
    { // composite Simpson's 1/3 rule for the first n-4 intervals plus Simpson's 3/8 rule for the last 3 intervals
        return simpson(n-3, f, dx) + simpson(4, &f[n-4], dx);
    }
}

double Integral::simpson(const int n, const double* const f, const double* const h)
{   
    // Simpson's rule for irregularly-spaced grid
    // The treatment for even number of grid points is the same as that of the regularly-spaced grid case.
    assert( n >= 2 );
    assert( std::all_of(h, h+(n-1), [](double x){return x > 0.0;}) );

    if (n == 4)
    {
        double w = h[0] + h[1] + h[2];
        return w / 12.0 * ( 2.0 + ((h[1]+h[2])/h[0]-1.0) * (h[2]/(h[0]+h[1])-1.0) ) * f[0]
            + std::pow(w,3) / 12.0 * (h[0]+h[1]-h[2]) / (h[0]*h[1]*(h[1]+h[2])) * f[1]
            + std::pow(w,3) / 12.0 * (h[2]+h[1]-h[0]) / (h[2]*h[1]*(h[1]+h[0])) * f[2]
            + w / 12.0 * ( 2.0 + ((h[1]+h[0])/h[2]-1.0) * (h[0]/(h[2]+h[1])-1.0) ) * f[3];
    }

    if (n == 2)
    {
        return 0.5 * h[0] * (f[0] + f[1]);
    }

    if (n % 2 == 1)
    {
        double sum = 0.0;
        for (int i = 0; i < n/2; ++i)
        {
            double hrp = h[2*i+1] / h[2*i];
            double hrm = h[2*i] / h[2*i+1];
            sum += (h[2*i+1] + h[2*i]) / 6.0 * ( (2.0-hrp)*f[2*i] + (2.0+hrp+hrm)*f[2*i+1] + (2.0-hrm) * f[2*i+2]);
        }
        return sum;
    }
    else
    {
        return simpson(n-3, f, h) + simpson(4, &f[n-4], &h[n-4]);
    }
}

void Integral::Gauss_Legendre_grid_and_weight(const int n, double *x, double *weights)
{
    assert( n >= 1 );

    double z = 0.0;
    double z1 = 0.0;
    double p1 = 1.0;
    double p2 = 0.0;
    double p3 = 0.0;
    double pp = 0.0;

    int half_grid_num = static_cast<int>((n+1)/2);
    for(int i = 1; i <= half_grid_num; i++)
    {
        z = cos(ModuleBase::PI * (i - 0.25) / (n + 0.5));

        while(true)
        {
            p1 = 1.0;
            p2 = 0.0;

            for(int j = 1; j <= n; j++)
            {
                p3 = p2;
                p2 = p1;
                p1 = ((2.0 * j - 1.0) * z * p2 - (j - 1.0) * p3) / j;
            }

            pp = n * (p2 - z * p1) / (1.0 - z*z);
            z1 = z;
            z = z1 - p1 / pp;

            if (std::abs(z - z1) < 1e-13) break;
        }

        x[i-1] = -z;
        x[n-i] = z;
        weights[i-1] = 2.0 / ((1.0 - z * z) * pp * pp);
        weights[n-i] = weights[i-1];
    }

}

void Integral::Gauss_Legendre_grid_and_weight(const double xmin, const double xmax, const int n, double *x, double *weights)
{
    assert( n >= 1 );

    double xl = (xmax - xmin) * 0.5;
    double xmean = (xmax + xmin) * 0.5;

    double z = 0.0;
    double z1 = 0.0;
    double p1 = 1.0;
    double p2 = 0.0;
    double p3 = 0.0;
    double pp = 0.0;

    int half_grid_num = static_cast<int>((n+1)/2);
    for(int i = 1; i <= half_grid_num; i++)
    {
        z = cos(ModuleBase::PI * (i - 0.25) / (n + 0.5));

        while(true)
        {
            p1 = 1.0;
            p2 = 0.0;

            for(int j = 1; j <= n; j++)
            {
                p3 = p2;
                p2 = p1;
                p1 = ((2.0 * j - 1.0) * z * p2 - (j - 1.0) * p3) / j;
            }

            pp = n * (p2 - z * p1) / (1.0 - z*z);
            z1 = z;
            z = z1 - p1 / pp;

            if (std::abs(z - z1) < 1e-13) break;
        }

        x[i-1] = xmean - xl * z;
        x[n-i] = xmean + xl * z;
        weights[i-1] = 2.0 * xl / ((1.0 - z * z) * pp * pp);
        weights[n-i] = weights[i-1];
    }

}

const double Integral::Lebedev_Laikov_grid110_x[110] = {
   1.000000000000000,  -1.000000000000000,   0.000000000000000,   0.000000000000000,   0.000000000000000,   0.000000000000000,   0.577350269189626,  -0.577350269189626,
   0.577350269189626,  -0.577350269189626,   0.577350269189626,  -0.577350269189626,   0.577350269189626,  -0.577350269189626,   0.185115635344736,  -0.185115635344736,
   0.185115635344736,  -0.185115635344736,   0.185115635344736,  -0.185115635344736,   0.185115635344736,  -0.185115635344736,   0.185115635344736,  -0.185115635344736,
   0.185115635344736,  -0.185115635344736,   0.185115635344736,  -0.185115635344736,   0.185115635344736,  -0.185115635344736,   0.965124035086594,  -0.965124035086594,
   0.965124035086594,  -0.965124035086594,   0.965124035086594,  -0.965124035086594,   0.965124035086594,  -0.965124035086594,   0.690421048382292,  -0.690421048382292,
   0.690421048382292,  -0.690421048382292,   0.690421048382292,  -0.690421048382292,   0.690421048382292,  -0.690421048382292,   0.690421048382292,  -0.690421048382292,
   0.690421048382292,  -0.690421048382292,   0.690421048382292,  -0.690421048382292,   0.690421048382292,  -0.690421048382292,   0.215957291845848,  -0.215957291845848,
   0.215957291845848,  -0.215957291845848,   0.215957291845848,  -0.215957291845848,   0.215957291845848,  -0.215957291845848,   0.395689473055942,  -0.395689473055942,
   0.395689473055942,  -0.395689473055942,   0.395689473055942,  -0.395689473055942,   0.395689473055942,  -0.395689473055942,   0.395689473055942,  -0.395689473055942,
   0.395689473055942,  -0.395689473055942,   0.395689473055942,  -0.395689473055942,   0.395689473055942,  -0.395689473055942,   0.828769981252592,  -0.828769981252592,
   0.828769981252592,  -0.828769981252592,   0.828769981252592,  -0.828769981252592,   0.828769981252592,  -0.828769981252592,   0.478369028812150,  -0.478369028812150,
   0.478369028812150,  -0.478369028812150,   0.878158910604066,  -0.878158910604066,   0.878158910604066,  -0.878158910604066,   0.478369028812150,  -0.478369028812150,
   0.478369028812150,  -0.478369028812150,   0.878158910604066,  -0.878158910604066,   0.878158910604066,  -0.878158910604066,   0.000000000000000,   0.000000000000000,
   0.000000000000000,   0.000000000000000,   0.000000000000000,   0.000000000000000,   0.000000000000000,   0.000000000000000
};

const double Integral::Lebedev_Laikov_grid110_y[110] = {
   0.000000000000000,   0.000000000000000,   1.000000000000000,  -1.000000000000000,   0.000000000000000,   0.000000000000000,   0.577350269189626,   0.577350269189626,
  -0.577350269189626,  -0.577350269189626,   0.577350269189626,   0.577350269189626,  -0.577350269189626,  -0.577350269189626,   0.185115635344736,   0.185115635344736,
  -0.185115635344736,  -0.185115635344736,   0.185115635344736,   0.185115635344736,  -0.185115635344736,  -0.185115635344736,   0.965124035086594,   0.965124035086594,
  -0.965124035086594,  -0.965124035086594,   0.965124035086594,   0.965124035086594,  -0.965124035086594,  -0.965124035086594,   0.185115635344736,   0.185115635344736,
  -0.185115635344736,  -0.185115635344736,   0.185115635344736,   0.185115635344736,  -0.185115635344736,  -0.185115635344736,   0.690421048382292,   0.690421048382292,
  -0.690421048382292,  -0.690421048382292,   0.690421048382292,   0.690421048382292,  -0.690421048382292,  -0.690421048382292,   0.215957291845848,   0.215957291845848,
  -0.215957291845848,  -0.215957291845848,   0.215957291845848,   0.215957291845848,  -0.215957291845848,  -0.215957291845848,   0.690421048382292,   0.690421048382292,
  -0.690421048382292,  -0.690421048382292,   0.690421048382292,   0.690421048382292,  -0.690421048382292,  -0.690421048382292,   0.395689473055942,   0.395689473055942,
  -0.395689473055942,  -0.395689473055942,   0.395689473055942,   0.395689473055942,  -0.395689473055942,  -0.395689473055942,   0.828769981252592,   0.828769981252592,
  -0.828769981252592,  -0.828769981252592,   0.828769981252592,   0.828769981252592,  -0.828769981252592,  -0.828769981252592,   0.395689473055942,   0.395689473055942,
  -0.395689473055942,  -0.395689473055942,   0.395689473055942,   0.395689473055942,  -0.395689473055942,  -0.395689473055942,   0.878158910604066,   0.878158910604066,
  -0.878158910604066,  -0.878158910604066,   0.478369028812150,   0.478369028812150,  -0.478369028812150,  -0.478369028812150,   0.000000000000000,   0.000000000000000,
   0.000000000000000,   0.000000000000000,   0.000000000000000,   0.000000000000000,   0.000000000000000,   0.000000000000000,   0.478369028812150,  -0.478369028812150,
   0.478369028812150,  -0.478369028812150,   0.878158910604066,  -0.878158910604066,   0.878158910604066,  -0.878158910604066
};

const double Integral::Lebedev_Laikov_grid110_z[110] = {
   0.000000000000000,   0.000000000000000,   0.000000000000000,   0.000000000000000,   1.000000000000000,  -1.000000000000000,   0.577350269189626,   0.577350269189626,
   0.577350269189626,   0.577350269189626,  -0.577350269189626,  -0.577350269189626,  -0.577350269189626,  -0.577350269189626,   0.965124035086594,   0.965124035086594,
   0.965124035086594,   0.965124035086594,  -0.965124035086594,  -0.965124035086594,  -0.965124035086594,  -0.965124035086594,   0.185115635344736,   0.185115635344736,
   0.185115635344736,   0.185115635344736,  -0.185115635344736,  -0.185115635344736,  -0.185115635344736,  -0.185115635344736,   0.185115635344736,   0.185115635344736,
   0.185115635344736,   0.185115635344736,  -0.185115635344736,  -0.185115635344736,  -0.185115635344736,  -0.185115635344736,   0.215957291845848,   0.215957291845848,
   0.215957291845848,   0.215957291845848,  -0.215957291845848,  -0.215957291845848,  -0.215957291845848,  -0.215957291845848,   0.690421048382292,   0.690421048382292,
   0.690421048382292,   0.690421048382292,  -0.690421048382292,  -0.690421048382292,  -0.690421048382292,  -0.690421048382292,   0.690421048382292,   0.690421048382292,
   0.690421048382292,   0.690421048382292,  -0.690421048382292,  -0.690421048382292,  -0.690421048382292,  -0.690421048382292,   0.828769981252592,   0.828769981252592,
   0.828769981252592,   0.828769981252592,  -0.828769981252592,  -0.828769981252592,  -0.828769981252592,  -0.828769981252592,   0.395689473055942,   0.395689473055942,
   0.395689473055942,   0.395689473055942,  -0.395689473055942,  -0.395689473055942,  -0.395689473055942,  -0.395689473055942,   0.395689473055942,   0.395689473055942,
   0.395689473055942,   0.395689473055942,  -0.395689473055942,  -0.395689473055942,  -0.395689473055942,  -0.395689473055942,   0.000000000000000,   0.000000000000000,
   0.000000000000000,   0.000000000000000,   0.000000000000000,   0.000000000000000,   0.000000000000000,   0.000000000000000,   0.878158910604066,   0.878158910604066,
  -0.878158910604066,  -0.878158910604066,   0.478369028812150,   0.478369028812150,  -0.478369028812150,  -0.478369028812150,   0.878158910604066,   0.878158910604066,
  -0.878158910604066,  -0.878158910604066,   0.478369028812150,   0.478369028812150,  -0.478369028812150,  -0.478369028812150
};

const double Integral::Lebedev_Laikov_grid110_w[110] = {
   0.048107465851397,   0.048107465851397,   0.048107465851397,   0.048107465851397,   0.048107465851397,   0.048107465851397,   0.123071735281670,   0.123071735281670,
   0.123071735281670,   0.123071735281670,   0.123071735281670,   0.123071735281670,   0.123071735281670,   0.123071735281670,   0.103191734088330,   0.103191734088330,
   0.103191734088330,   0.103191734088330,   0.103191734088330,   0.103191734088330,   0.103191734088330,   0.103191734088330,   0.103191734088330,   0.103191734088330,
   0.103191734088330,   0.103191734088330,   0.103191734088330,   0.103191734088330,   0.103191734088330,   0.103191734088330,   0.103191734088330,   0.103191734088330,
   0.103191734088330,   0.103191734088330,   0.103191734088330,   0.103191734088330,   0.103191734088330,   0.103191734088330,   0.124945096872513,   0.124945096872513,
   0.124945096872513,   0.124945096872513,   0.124945096872513,   0.124945096872513,   0.124945096872513,   0.124945096872513,   0.124945096872513,   0.124945096872513,
   0.124945096872513,   0.124945096872513,   0.124945096872513,   0.124945096872513,   0.124945096872513,   0.124945096872513,   0.124945096872513,   0.124945096872513,
   0.124945096872513,   0.124945096872513,   0.124945096872513,   0.124945096872513,   0.124945096872513,   0.124945096872513,   0.120580249028528,   0.120580249028528,
   0.120580249028528,   0.120580249028528,   0.120580249028528,   0.120580249028528,   0.120580249028528,   0.120580249028528,   0.120580249028528,   0.120580249028528,
   0.120580249028528,   0.120580249028528,   0.120580249028528,   0.120580249028528,   0.120580249028528,   0.120580249028528,   0.120580249028528,   0.120580249028528,
   0.120580249028528,   0.120580249028528,   0.120580249028528,   0.120580249028528,   0.120580249028528,   0.120580249028528,   0.121830917385521,   0.121830917385521,
   0.121830917385521,   0.121830917385521,   0.121830917385521,   0.121830917385521,   0.121830917385521,   0.121830917385521,   0.121830917385521,   0.121830917385521,
   0.121830917385521,   0.121830917385521,   0.121830917385521,   0.121830917385521,   0.121830917385521,   0.121830917385521,   0.121830917385521,   0.121830917385521,
   0.121830917385521,   0.121830917385521,   0.121830917385521,   0.121830917385521,   0.121830917385521,   0.121830917385521
};
}
